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Abstract 

Second order Sobolev metrics on the space of regular unparametrized planar curves 
have several desirable completeness properties not present in lower order metrics, but nu¬ 
merics are still largely missing. In this paper, we present algorithms to numerically solve 
the initial and boundary value problems for geodesics. The combination of these algo¬ 
rithms allows to compute Karcher means in a Riemannian gradient-based optimization 
scheme. Our framework has the advantage that the constants determining the weights 
of the zero, first, and second order terms of the metric can be chosen freely. Moreover, 
due to its generality, it could be applied to more general spaces of mapping. We demon¬ 
strate the effectiveness of our approach by analyzing a collection of shapes representing 
physical objects. 


Introduction 


Unparametrized curves arise naturally in shape analysis and its many applications, including 
medical imaging [El, E3, El], object tracking [ED, EID], computer animation [Q, O], computer 
aided design [D], speech recognition [ED], analysis of bird migration patterns and hurricane 
paths [E3], biology [0,113], and many other fields [□, D]. In many instances, the rationale 
for identifying curves differing only by a reparametrization is that the curves represent the 
boundaries of physical shapes. Shapes can be analyzed mathematically by endowing the 
space of shapes with a Riemannian metric. Riemannian shape analysis has developed into 
an active field of research by now. 

To do statistics on shape space, robust and efficient implementations of the boundary and 
initial value problems for geodesics are needed. Unfortunately, the arguably simplest metric 
on shape space is degenerate: it is well-known that the Riemannian distance induced by the 
L^-metric vanishes on spaces of parametrized and unparametrized curves [ED]. The discovery 
of this degeneracy led to an investigation of Sobolev metrics of higher order [ED, ED, E3]. 
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For the important class of first order metrics numerics are well developed by now, but 
a restriction on the parameters of the metric must be imposed [O, 113]. The situation is 
drastically different for second order metrics, despite the fact that they enjoy better com¬ 
pleteness properties [H, 0]. On the positive side, the geodesic boundary value problem under 
second order Finsler metrics on the space of BV^ curves was implemented numerically in 
[IZ3]. Moreover, on the space of parametrized curves, there are numerics for the boundary 
and initial value problems for geodesics under second order Sobolev metrics [i, □]. On 
spaces of unparametrized curves numerics for second order Sobolev metrics are, however, 
still lacking. This is an important challenge and the topic of this paper. 

We present a numerical implementation of the initial and boundary value problem for 
geodesics of planar, unparametrized curves under second order Sobolev metrics. Our imple¬ 
mentation is based on a discretization of the Riemannian energy functional using B-splines. 
The boundary value problem for geodesics is solved by a gradient descent on the set of dis¬ 
cretized paths and the initial value problem by discrete geodesic calculus [E3]. Our approach 
is general in that it involves no restriction on the parameters of the metric and allows to factor 
out rigid transformations. 

We tested our implementation on a dataset of shapes representing various groups of simi¬ 
lar physical objects [113]. We were able to reconstruct the groups using agglomerative cluster¬ 
ing based on the matrix of pairwise geodesic distances between the shapes. We also studied 
within-group shape variations by performing nonlinear principal component analysis using 
the Karcher means of the groups as base points. We visualized the results by solving the 
geodesic equation forward and backward in time along the principal directions. 

Our algorithms performed well on the shapes in the database. Some care was needed 
to avoid singularities in the spline representations of the contours extracted from the binary 
images and of the paths initializing the gradient descent of the energy functional. 

In future work, our framework could be applied to other spaces of mappings like manifold¬ 
valued curves, embedded surfaces, or more general spaces of immersions (see [□, □] for 
details and [Q] for a general overview). A rigorous convergence analysis of the proposed 
discretization remains to be done. 


2 Mathematical background 

We extend the exposition of [□] to unparametrized curves. The space of smooth, regular, 
parametrized curves is defined as 

Imm(5',K'') = {ceC°°(5SM''):V0e5‘,C0(0) +0} , 

where Imm stands for immersion. Imm(5',K'^) is an open subset of the Frechet space 
and as such itself a Frechet manifold. The tangent space 7’cImm(5',K'^) at any 
curve c is the vector space Geometrically, tangent vectors at c are K^-valued 

vector fields along c. 

We denote the Euclidean inner product on by (•,•). Moreover, for any fixed curve 
c, we denote differentiation and integration with respect to arc length by Dg = and 

ds = \cg\d9, respectively. 

Definition 2.1. Second order Sobolev metrics on Imm(5^,K'^) are Riemannian metrics of 
the form 

Gc{h,k) = J ^ aQ{h,k) + a\{Dsh,Dsk) +a2{D^h,D^k)ds , 
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where ao,a 2 > 0, ai > 0 are constants and h,ke tangent vectors. 

Remark 2.2 (Choice of the constants). In the definition of Sobolev metrics, there is freedom 
in the choice of the relative weighting of the L^-, H^- and //^-parts. This choice influences 
the geometry of the space of curves and should, in each application, be informed by the data 
at hand. 

The reparametrization group is the diffeomorphism group of the circle, 

Diff(5') = {(p€C“(5',5^) : ^bijective}, 

which is an infinite-dimensional regular Frechet Lie group [EB]. Reparametrizations act on 
curves by composition from the right, {c,(p) co (p. The shape space of unparametrized 
curves is the orbit space 5(5'= Imm(5',K^)/Diff(5*) of this group action.' 

Theorem 2.3 ([EU], Sect. 1.5). The space 5(5',M'^) is a Frechet manifold and the base 
space of the principal fibre bundle 

;r:Imm(5',K'') ^5(5',K"'), ccoDiff(5'), 

with structure group Diff(5'). A Sobolev metric G on Imm(5',M‘^) induces a metric on 
5(5', such that the projection K is a Riemannian submersion. 

Remark 2.4 (Invariance of the metric). Thm. 2.3 hinges on the invariance of Sobolev met¬ 
rics with respect to reparametrizations. Sobolev metrics are also invariant with respect to 
translations and rotations, but in general not to scalings. The lack of scale invariance can be 
addressed by introducing weights depending on the length 4 of the curve c as follows: 

Gcihfi) = ‘^^{h,k) + ‘^{D,h,D,k) + a2tc{D]h,D]k)ds. 

Deformations of curves are smooth paths c: [0,1] ^ Imm(5' ,]R‘^). Their velocity is q, 
the subscript t denoting differentiation. The length of a path c is 

L(c) = ^ ^JG,(t){ct{t),c,{t))dt. 

The distance between two shapes K{cf) and 7r(ci) in 5(5',with respect to the metric 
induced by G is the infimum over all paths between co and the orbit ci oDiff(5'), i.e., 

dist(co,ci)= inf L{c). 

c(0)=Co 

e(l)6e|oDiff(s‘) 

Geodesics on Imm(5',IR'^) are critical points of the energy functional 

1 r ^ 

= G,(^t^{c,{t),c,{t))dt. 

The spaces 5(5' and Imm(5' ,K‘^) are related by a Riemannian submersion and geodesics 
on 5(5',K'^) can be lifted to horizontal geodesics on Imm(5',K‘^). Conversely, the projec¬ 
tion of a length-minimizing path between cq and ci oDiff(5' is a geodesic on 5(5', K"^). 

The space 5(5',M‘^) equipped with a Sobolev metric possesses some nice completeness 
properties, which are summarized in the following theorem. 


'More precisely, we define iS(S*,R'') = Immf(S*,R'')/Diff(S*), where Immf(5’,R'') is the subset of free 
immersions, i.e., those upon which Diff(S') acts freely. This is an open and dense subset of Imm(S*,R‘'). While 
important for theoretical reasons, this restriction has no influence on the practical applications of Sobolev metrics. 
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Theorem 2.5 (Minimizing geodesics, [H, 0]). Let ao,a 2 > 0 and ai>0. Then, given two 
curves co,ci in the same connected component there exists a minimizing 

geodesic connecting them. Furthermore, there exists a minimizing geodesic connecting the 
shapes 7r(co) and ?r(ci) in 5(5*,M'^). 

Remark 2.6 (Elastic metrics). Closely related to the Sobolev metrics described here is the 
family of elastic metrics [IZH], which in the planar case is given by 

Gc(h,k) = f a^{Dsh,n){Dsk,n) +b^{Dsh,v){Dsk,v)ds . 

Js^ 

Here are constants and v,n denote the unit tangent and normal vectors to c. Two special 
cases deserve to be highlighted; for a = 1, = 2 [Q] and a = b [E3] there exist nonlinear 

transforms, the square root velocity transform and the basic mapping, that greatly simplify 
the numerical computation of geodesics. Both of these metrics have been applied to a variety 
of problems in shape analysis. We note that the elastic metric with a = b corresponds to a first 
order Sobolev metric as in Def. 2.1 with ao = a 2 = 0 and ai = a^ = b^. As it has no L^-part, it 
is a Riemannian metric only on the space of curves modulo translations. 

3 Numerical implementation 

3.1 Discretization 

We discretize curves using B-splines; c= djCj{0), where Cj are the B-splines of degree 
ng, defined on a uniform periodic knot sequence, with all knots of multiplicity one. Observe 
that Cj € C"®“*([0,27r]). A path of curves can then be represented using tensor product 
B-splines, i.e., 

N, Ng 

c{t,e) = Y,t.di.jB,{t)Cj{e). (1) 

;=!>! 

Here Bi are the B-splines defined on the interval [0,1], with uniform knots and full multi¬ 
plicity at the end points, Bj e C"'“*([0,1]). This implies that the boundary curves are given 
by c(0, 9 ) = dijCj(e) and c( 1,0) = dN,jCj(e). 

Under the identification of 5* with IR/[0,27r], diffeomorphisms 5* can be written 

as Xj/ = id+0, where (j) is a periodic function. We choose to discretize (j){9) = (j)iDi{9), 

where D, are B-splines of degree n^, defined on a uniform periodic knot sequence, similarly 
to Cj. The identity can be written in a B-spline basis using the Greville abscissas i.e., 

id = ^iDj. Due to the positivity of B-splines, the condition Xj/' > 0 ensuring that i// is a 
diffeomorphism takes the form 

( 2 ) 

To speed up convergence, we introduce an additional variable a € K representing constant 
shifts of the reparametrization. The resulting redundancy is eliminated by the constraint 

T.^i = 0. ( 3 ) 

(-1 

To compute the energy of paths (1), we use Gaussian quadrature on each interval between 
subsequent knots. Evaluating paths (and their derivatives) at the quadrature points is a simple 
multiplication of the spline collocation matrix with the vector of control points. 
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Figure 1; Selection of shapes from the dataset [D3]. 


3.2 Geodesics and Karcher means 

From now on we work with plane curves (d = 2). The boundary value problem for geodesics 
consists of minimizing the discretized energy (2) over all paths c: [0,1] x [0,27r] ^ map¬ 
pings \j/ = id+0:5^ ^ 5*, shifts of the reparametrization a 6 K, and rotations around the 
origin by an angle j3 e [0,2;r), subject to the constraints (2), (3), and 

c(0,-) = co(-), c(l,-) =/?jg(ci (!//(•)-«)+0, 

where co,ci: [0,2;r] ^ are given boundary curves. This is a hnite dimensional constrained 
optimization problem, which we solve using Matlab’s interior point method fmincon. We 
achieved major performance improvements by fine-tuning the implementations of the gradi¬ 
ent and hessian of the energy functional. 

In the initial value problem for geodesics, the initial value and initial velocity of the 
sought geodesic are given. As described in Sect. 2, geodesics on can be lifted 

to horizontal geodesics on Imm(5',K‘^). Thus, solving the geodesic initial value problem 
for unparametrized curves reduces to solving the problem for parametrized curves with hor¬ 
izontal initial velocities. To solve the latter problem, we use the time-discrete variational 
geodesic calculus of [IZ3] as described in [□]. 

The Karcher mean c of a set {ci,..., c„ } of curves is dehned as the minimizer of 

";=i 

It can be computed by iteratively solving initial and boundary value problems for geodesics. 
We refer to [E3] for more details. 


4 Numerical examples 

4.1 Data acquisition and setup 

We tested our implementation on a dataset of shapes collected by the Computer Vision Group 
at Brown university [O]. The dataset consists of black and white images of physical objects. 
It is natural to represent the objects by their boundaries using unparametrized curves. In 
addition to factoring out reparametrizations, we also factor out translations and rotations 
because we are not interested in the position of the curves in space. Some of the resulting 
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Figure 2: Geodesics between a fish and a tool in the space of unparametrized curves. The 
metric parameter ai is increased by a factor 10 in the second, a factor 100 in the third, and a 
factor 1000 in the fourth column. The corresponding geodesic distances are 138.02, 162.55, 
246.46 and 468.74. Note that since we also optimize over translations and rotations of the 
target curve, the position in space varies. 


curves are depicted in Fig. 1 . We used splines of degree ng =n^ = 3 with Ng = 60 and = 20 
controls in space and of degree n, = 2 with Nt = 20 controls in time. 

To solve the boundary value for geodesics as described in Sect. 3.2, one has to construct 
an initial homotopy between the given boundary curves. If the curves differ only by a small 
deformation, this is unproblematic because the linear path connecting the curves can be used. 
For larger deformations, however, the linear path often has self-intersections, changing the 
winding number, which can lead to convergence problems during the energy minimization. 
Our solution was to construct homotopies by deforming the initial curve to a circle and the 
circle to the target curve using linear paths. This worked well for all examples presented 
in this article. It also led to the same optima in all cases where the linear path was free of 
self-intersections. In future work, we plan to investigate further ways of creating good and 
natural initial paths. 

The choice of parameters ag, ai, and 02 of the Riemannian metric can have a large 
infiuence on the resulting optimal deformations, as can be seen in Fig. 2. In this article, 
we used the following ad-hoc strategy for choosing the constants: we computed the average 
L^-, and //^-contributions £ 0 , £ 1 , £2 to the energy of linear paths between each pair of 
curves in the dataset. Then we normalized ag to 1 and chose ai and a 2 such that 

agEg : ai£i : aiEi =1:1:1 and £ = agEg + a\E\ + a2£2 = 100. 

This resulted in the parameter values a\ = 250 and a 2 = 0.004. It remains open how to best 
choose the constants depending on the data under consideration. 

4.2 Comparison to non-elastic metrics 

Riemannian metrics on spaces of unparametrized curves are often called elastic metrics, 
as they allow both for bending and stretching of the curve. In the elastic case, solving 
the boundary value problem for geodesics involves optimizing over the Diff(5*)-orbit of 
the initial or final shape. This is a computationally expensive and difficult task since the 
diffeomorphism group is infinite-dimensional. 

An alternative and simpler approach is to parametrize the curves by unit-speed and to 
calculate geodesics in the space of parametrized curves. This could be called a non-elastic 
approach. (Of course, one might wish to also factor out rigid transformations and constant 
shifts of the parametrization, but this is much simpler because these groups are finite dimen¬ 
sional.) 
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Figure 3: Geodesics in the space of unparametrized (1st, 3rd) versus parametrized (2nd, 4th) 
curves modulo rotations and translations. Note that since we also optimize over translations 
and rotations of the target curve, the curves in the first two panels are aligned differently. 
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Figure 4: The matrix of geodesic distances between the shapes in Fig. 1, visualized using 
multi-dimensional scaling in two and three dimensions. The labels are: fish O, sting rays ☆, 
bunnies tools +. 


Our experiments suggest that using the more involved elastic approach pays off. The 
two approaches yield different results, and in particular in cases involving large amounts of 
stretching the geodesics found using the first approach appear more natural. However, in 
cases involving mainly bending of the curves the results are very similar (see Fig. 3). 


4.3 Clustering and principal component analysis 

We investigated if pairwise geodesic distances can be used to cluster shapes into meaningful 
groups. To this aim, we calculated all pairwise distances between the 33 shapes presented 
in Fig. 1. Solving the corresponding 528 boundary value problems took about two hours 
on a 3GHz processor with four cores. The resulting distance matrix is visualized in Fig. 4 
using multi-dimensional scaling; the plot suggests that objects of the same group lie close 
together. Indeed, agglomerative clustering with 4 clusters reproduces exactly the subgroups 
of the database. 

Finally, we studied within-group variations using non-linear principal component anal¬ 
ysis. To this aim, we first computed the Karcher mean of each group. The corresponding 
optimization problem (4) was solved using a conjugate gradient method, as implemented 
in the Manopt library [0], on the finite-dimensional spline approximation of the Riemannian 
manifold of curves. The mean shapes for the groups of fish and humans can be seen in Fig. 5. 
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Figure 5: First column: Karcher means (bold) of the groups of fish and humans. Second and 
third column: geodesics from the mean in the first and second principal direction at the times 
-3,-2,... ,2,3; the bold curve is the mean. 


Next, we represented each shape in the group by the initial velocity from the mean c using 
the inverse of the Riemannian exponential map. We then performed a principal component 
analysis with respect to the inner product Gj on the set of initial velocities. In the group 
of human figures, the first three eigenvalues capture 67%, 22%, and 6% of within-group 
variation. In the group of hsh, the first three eigenvalues capture only 40%, 25%, and 16% 
of within-group variation. Geodesics from the mean in the directions of the hrst two prin¬ 
cipal directions can be seen in Fig. 5. In the group of humans the first principal direction 
encodes bending of the arms and legs, whereas the second direction reflects stretching in the 
extremities. 


5 Conclusions 

In this article we developed a numerical framework for solving the initial and boundary value 
problem for geodesics of planar, unparametrized curves under second order Sobolev metrics. 
We tested our implementation on a dataset of shapes representing various groups of similar 
physical objects and obtained good experimental results. In future work, we plan to apply our 
algorithms to datasets of real-world medical data, prove rigorous convergence results for the 
discretization, and extend the framework to other spaces of mappings like manifold-valued 
curves and embedded surfaces. 
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